#Search-Tissue

#library packages
library(ggplot2)
library(dplyr)
library(data.table)

library(optparse)
library(stringr)

option_list <- list(
  make_option("--i", default = "", type = "character", help = "input file"),
  make_option("--m", default = "", type = "character", help = "meta file"),
  make_option("--c", default = "", type = "character", help = "color file"),
  make_option("--g", default = "", type = "character", help = "group"),
  make_option("--p", default = "", type = "character", help = "phenotype names use for sample filter"),
  make_option("--d", default = "", type = "character", help = "dataset names use for sample filter"),
  make_option("--model", default = "", type = "character", help = "model names use for sample filter"),
  make_option("--l", default = "", type = "character", help = "log"),
  make_option("--o", default = "", type = "character", help = "output png file")
)
opt <- parse_args(OptionParser(option_list = option_list))


#import source data
expr <- fread(opt$i) %>%
  as.data.frame()
rownames(expr) <- expr$gene
gene <- expr$gene
expr <- expr[, -1]
meta <- read.table(opt$m, header = T, stringsAsFactors = F, sep = '\t')
mypal <- read.table(opt$c)[, 1] %>%
  as.character()

#change strings into factors #修改
meta$Phenotype <- factor(meta$Phenotype, levels = c('Fetal', 'Normal', 'Fatty Liver', 'NAFLD', 'Viral Hepatitis', 'Fibrosis',
                                                    'Control', 'HCC', 'HB', 'Ductular Proliferation', 'IPNB', 'ICC'))
meta <- arrange(meta, Phenotype, Model)
meta$Model <- factor(meta$Model, levels = c('Fetal', 'Normal', 'Control', unique(meta$Model)[4:length(unique(meta$Model))]))
meta <- arrange(meta, Model, Duration)
meta$Model_duration <- factor(meta$Model_duration, levels = unique(meta$Model_duration))

#dim parameters
group <- opt$g #choose one from ('Phenotype','Model','Model_duration')
phenotype <- opt$p %>%
  str_split(";") %>%
  unlist()
dataset <- opt$d %>%
  str_split(";") %>%
  unlist()
opt$model
model <- opt$model %>%
  str_split("@@@") %>%
  unlist()
model
filter <- list(phenotype = phenotype, dataset = dataset, model = model)
logscale <- opt$l

#extract data
meta.filter <- meta[(meta$Phenotype %in% filter$phenotype) &
                      (meta$Dataset %in% filter$dataset) &
                      (meta$Model_show %in% filter$model),]

tmpFilter<- meta[(meta$Phenotype %in% filter$phenotype) &
                   (meta$Model_show %in% filter$model),]
filter$model
nrow(tmpFilter)
meta.filter

expr.gene <- expr[gene, meta.filter$Run]
meta.filter <- mutate(meta.filter, expr = as.numeric(expr.gene), log2expr = log2(as.numeric(expr.gene) + 1),
                      group = meta.filter[, which(colnames(meta.filter) == group)])

#calculate freq
freq <- as.data.frame(table(meta.filter$group))
freq <- freq[freq$Freq != 0,]
if (group == 'Phenotype') {
  meta.filter$group <- factor(meta.filter$group, levels = freq$Var1, labels = paste(gsub(' ', '\n', freq$Var1), '\n(', freq$Freq, ')', sep = ''))
} else if (group %in% c('Model', 'Model_duration')) {
  meta.filter$group <- factor(meta.filter$group, levels = freq$Var1, labels = paste(freq$Var1, ' (', freq$Freq, ')', sep = ''))
}

#plot #修改
if (group %in% c('Phenotype', 'Model')) {
  if (logscale == 'Yes') {
    p <- ggplot(meta.filter) +
      geom_boxplot(aes(x = group, y = log2expr, color = group), outlier.size = 0.8) +
      scale_color_manual(values = mypal) +
      theme_bw() +
      theme(legend.position = 'none') +
      labs(y = 'Log2(TPM+1)', x = group, title = gene) +
      theme(axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1, size = 7), axis.title.x = element_blank(),
            plot.title = element_text(hjust = 0.5))
  } else {
    p <- ggplot(meta.filter) +
      geom_boxplot(aes(x = group, y = expr, color = group), outlier.size = 0.8) +
      scale_color_manual(values = mypal) +
      theme_bw() +
      theme(legend.position = 'none') +
      labs(y = 'TPM', x = group, title = gene) +
      theme(axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1, size = 7), axis.title.x = element_blank(),
            plot.title = element_text(hjust = 0.5))
  }
} else if (group == 'Model_duration') {
  if (logscale == 'Yes') {
    p <- ggplot(meta.filter) +
      geom_boxplot(aes(x = group, y = log2expr, color = Model), outlier.size = 0.8) +
      scale_color_manual(values = mypal) +
      theme_bw() +
      theme(legend.position = 'none') +
      labs(y = 'Log2(TPM+1)', x = group, title = gene) +
      theme(axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1, size = 7), axis.title.x = element_blank(),
            plot.title = element_text(hjust = 0.5))
  } else {
    p <- ggplot(meta.filter) +
      geom_boxplot(aes(x = group, y = expr, color = Model), outlier.size = 0.8) +
      scale_color_manual(values = mypal) +
      theme_bw() +
      theme(legend.position = 'none') +
      labs(y = 'TPM', x = group, title = gene) +
      theme(axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1, size = 7), axis.title.x = element_blank(),
            plot.title = element_text(hjust = 0.5))
  }
}
png(file = opt$o, width = 2000, height = 1200, res = 300)
print(p)
dev.off()

#export csv
meta.filter$Size <- str_extract(as.character(meta.filter$group), '\\d+')
meta.filter1 <- select(meta.filter, Run, Phenotype, Dataset, Model_duration, Model, Size, expr, log2expr)
write.csv(meta.filter1, file = paste('gene.csv', sep = ''), row.names = F)
